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Abstract 

We use molecular dynamics simulations to understand the mechanisms of polarization switching 
in ferroelectric BaTi03 achieved with external electric field. For tetragonal and orthorhombic fer- 
roelectric phases, we determine the switching paths, and show that polarization rotation through 
intermediate monoclinic phases (a) facilitates switching at low fields (b) is responsible for a sharp 
anisotropy in polarization switching. We develop understanding of this through determination of 
detailed electric field-temperature phase diagrams, that are fundamental to technological applica- 
tions based on electromechanical and switching response of ferroelectrics. 

PACS numbers: 
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Ferroelectric (FE) materials exhibit spontaneous electric polarization whose magnitude 
and direction depend sensitively on temperature, pressure and electric field 1|. They are 
technologically very important because of two of their properties: (a) large electromechani- 
cal coupling that is exploited in their applications as sensors and actuators in micro-electro- 
mechanical systems [l| , and (b) switchability of their polarization from one state to another 
with applied field, that makes them useful in non-volatile memory devices like ferroelectric 
random access memories (FeRAMs) . Fundamental understanding of the mechanisms re- 
sponsible for these properties in the bulk crystals is essential to development of ferroelectrics 
with improved properties and their use in nano-scale devices. 

Direction of spontaneous polarization is typically along a crystallographic direction. 
For example, tetragonal, orthorhombic and rhombohedral ferroelectric states are charac- 
terized by polarization along [001], [110] and [111] directions, respectively. Polarization 
rotation 3) from a rhombohedral state towards a tetragonal one through an intermediate 
monoclinic phase was shown to be responsible for an ultrahigh piezoelectric coupling ob- 
served experimentally^. Monoclinic phases, characterized by the polarization vector (P) 
along low symmetry directions, were carefully analyzed and shown to be relevant to the 
large piezoelectricity in 92%PbZn 1/ / 3 Nb2/30 3 -8%PbTi03 js] and also found in another tech- 
nologically important ferroelectric, Pb(Zr 3; Tii_ 3 ;)03[6j. 

Interestingly, piezoelectric response of BaTiOa (a simple, classic ferroelectric) was shown 
to be enhanced significantly by crystallographically engineering single crystals through ap- 
plied DC biasM, ^1. This involves nonlinear response of BaTiOs and is linked with electric 
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field induced structural phase transitions [4], [5|, 12[ in ferroelectrics. Polarization switching, 
crucial to memory applications, is another closely related nonlinear phenomenon whose mi- 
croscopic mechanism is yet to be uncovered. Understanding of these nonlinear phenomena 
within a single picture is facilitated by knowledge of electric field-temperature phase dia- 
grams of a ferroelectric. 

Using phenomenological Landau- Ginzburg-Devonshire theory, electric field-temperature 
(E — T) phase diagrams have been studied for single crystal BaTiO"3, showing evidence for 



monoclinic phases 



101 ] . Similarly, it was shown in general by Vanderbilt and Cohen 



ll| that 



the extension of the sixth order Devonshire theory to eighth and twelfth orders in polarization 
as the order parameter was necessary to explain the presence of stable monoclinic and 
triclinic phases, respectively. 
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Here, we present molecular dynamics simulations carried out for a comprehensive inves- 
tigation of polarization switching in BaTiC>3 through determination of (a) its temperature 
dependence, (b) paths followed during switching with field along different directions, and (c) 
detailed electric field-temperature phase diagrams. This permits identification of different 
monoclinic phases that are stabilized as a function of electric field and temperature, and are 
relevant to both ultra-high piezoelectric response and polarization switching properties. 

We use a first-principles effective Hamiltonian in classical molecular dynamics simulations 
that (1) capture most nonlinearities by construction, and (2) include thermal fluctuations in 
polarization and strain. A similar approach has been used earlier jl2j] in finding polarization 
paths as a function of applied field at a fixed temperature. An effective Hamiltonian H e Q is 
a Taylor series expansion of the energy surface around the high symmetry cubic perovskite 
structure, written in terms of the low energy local degrees of freedoms (DoFs), three acoustic 
phonon modes (which capture the physics of inhomogeneous strain) and three low energy 
(soft) optical phonon modes (which give local dipole moments). In order to include effects 
of an external electric field on the system, a term — Z* S-u{R) where Z* and u(R) are the 
Born effective charge and displacement associated with the soft-mode variable of the unit 

n r 

cell at R, is added to H eS . We use parameters in H e g as determined in Refs |14| and [15] 
from first-principles density functional theory calculations with local density approximation 
(LDA). Use of LDA causes an underestimation of the lattice parameter, hence a negative 
pressure of -5 GPa is applied to all systems simulated jl5|] here. 

Mixed space molecular dynamics [la, 113] (MD) is used here to determine finite temper- 
ature properties of H e g- with periodic boundary conditions. The FERAM code used in 
our simulations is described in detail in Ref 19| and the program can be downloaded from 



http://lotto.sourceforge.net/feram In order to simulate polarization hysteresis, the 



temperature is first increased or decreased in steps of ~30 K to achieve an equilibrium at 
the temperature of the hysteresis simulation. Then, at fixed temperature, the electric field 
is switched on from zero to a positive value E max , from which it is decreased to -E max in 
steps of about 10 kV/cm to 100 kV/cm. Subsequently, it is again increased finally to E max 
to obtain a complete hysteresis loop. At each temperature (and electric field), the system 
is thermalized for initial 70,000 timesteps and then averaging is performed over the next 
30,000 timesteps, amounting to a simulation period of 0.2 ns. 

For the high temperature cubic paraelectric phase, we simulate response to electric field 
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TABLE I: Parameters in the Langevin function used in fitting P(E) data of the paraelectric phase 
at different temperatures. 



T(K) 


fci (C/m 2 ) 


k 2 (CmxlO- 5 ) 


400 


0.12 


25.4 


370 


0.12 


41.5 


340 


0.11 


101.8 



applied along [001] direction (symmetry would give the same response to fields along other 
directions). The Langevin function 



L{E) = hicothffi - g] ( 

describes fairly well the P(E) curves of the paraelectric phase (see Fig la and fitting param- 
eters ki and k 2 in Table 1). While ki, which relates to the saturated polarization remains 
roughly constant, k 2 , which is proportional to the dielectric constant of the paraelectric 
phase, diverges as T approaches T C (T <-> C). Almost perfect fits to P(E) curves of the 
paraelectric phase with just two parameters and their expected temperature dependence 
validates our MD simulations to study further these nonlinear hysteresis properties of the 
polar phases. 

For polar phases, we simulate response to field applied (a) along the respective polar axis 
and (b) along other high symmetry crystallographic directions (eg., along [101] direction 
of the tetragonal phase where P points along [001]). Generally, the magnitude of field at 
which switching occurs, S c , reduces as the temperature increases upto the transition to 
higher temperature phase (see results for tetragonal and rhombohedral phases in Fig lb and 
1c). 



1) 



We now focus on the direction-dependence of polarization switching in rhombohedral phase. 
We find that the coercive field S c when applied along [001] direction is almost a factor of 5 
smaller than that when applied along [111] direction. While the reduction in S c with T can 
be readily understood through the concept of lowering of Landau free energy barrier near 
the transition, a strong dependence of £ c on its direction is related to the fact that transition 
state itself could be different. 
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FIG. 1: (Color online) Average value of polarization P 2 as a function of electric field E for (a) 
cubic paraelectric and (b) tetragonal ferroelectric phases when electric field is applied along [001] 
direction, (c) Coercive field 8 C as a function of temperature for fields along [111], [101] and [001] 
directions respectively, applied to rhombohedral phase, (d) Cohen's Cube depicting the high 
symmetry tetragonal (T), orthorhombic (O) and rhombohedral (R) phases and monoclinic M^4, 
Mb and Mc phases. 

In order to investigate the transition states relevant to switching, we examine the paths 
followed by the system during polarization switching. When E is parallel to P, the path is 
linear (P to -P) by symmetry. We use the following nomenclature to describe the different 
phases encountered along the polarization paths (Fig. 2): T, O and R denote the tetragonal, 
orthorhombic and rhombohedral phases respectively, and, Ma, Mb and Mc denote the three 



different monoclinic phases as defined in Ref. 



In the M c phase, P points along [Ouv] 



(u^v), and in the Ma and Mb phases P points along [uuv], with u<v and u>v, respectively. 
The monoclinic phases are also illustrated with the Cohen's cube[3j (see Fig. Id). 

With field along [001] and [111] directions, switching in the orthorhombic phase (at 240 K) 
occurs through monoclinic phases Mc and Mb that appear along the path at 25 kV/cm and 
15\/3 kV/cm, respectively. The transition from orthorhombic to the respective monoclinic 
phases is marked by (1) a jump in the polarization path from point 1 to 2 in Fig. 2a and 2b, 
and (2) a change in the slope of the curve of eigenvalues rj of strain tensor as a function of 
electric field (Fig. 3a and 3b). The second crossover (Mc— >T and Mb^-R) manifests as a 
change in the slope of rj at ~ 60 kV/cm (Fig. 3a) and ~ 50 kV/cm (Fig. 3b), respectively. 
The sequence of phases, O— >M C ^T (Fig. 2a) and O— >M B — >R (Fig. 2b) can be readily 
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FIG. 2: Paths of polarization traced during switching when field is applied along (a) [001] and (b) 
[111] directions to orthorhombic phase at 240 K and (c) [101] and (d) [111] directions to tetragonal 
phase at 290 K. The numbers mark the sequence of phases (nomenclature in the text) followed by 
P. 

rationalized through continuous paths along the Cohen's cube (see Fig. Id). 

Along the polarization paths for tetragonal phase (at 290 K), the monoclinic phases Mc 
and M A emerge when fields point along [101] (Fig. 2c) and [111] (Fig. 2d) directions, respec- 
tively. The transitions from tetragonal to the monoclinic phases are evident in the changing 
slopes at ~ 40 kV/cm of rj (Fig. 3c and 3d). Only at higher fields ~ 100 kV/cm, close to 
E max , the phases become orthorhombic and rhombohedral, respectively. The sequence of 
phases T— ^Mc^O and T— >M^^R also follows from Cohen's cube (Fig. Id). 

At a field of 25 kV/cm along [001] direction of the orthorhombic phase, there occurs a 
jump (from point 6 to 7 in Fig. 2a) along the path of switching which corresponds to a 
rotation from M c to M C - Similarly when the field is 15\/3 kV/cm along [111] direction, 
a jump from point 5 to 6 (see Fig 2b) corresponds to a rotation from O to O. Both 
the switching fields are lower compared to S c of 40^ kV/cm along [101] direction. In the 
tetragonal phase as well, switching from T to T occurs at the field of 25a/2 kV/cm along [101] 
direction, as revealed by jump from point 7 to 8 in Fig. 2c, and at the field of lb\f?> kV/cm 
along [111] direction as revealed by a jump from point 5 to 6 in Fig. 2d; both £ c 's are much 
lower than the coercive field of S c = 60 kV/cm along [001] direction. Thus, our finding for 
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FIG. 3: (Color online) Eigenvalues rj of strain tensor as a function of field applied along (a) [001] 
(b) [111] directions to the orthorhombic phase at 240 K and (c) [101] (d) [111] directions to the 
tetragonal phase at 290 K. 

the rhombohedral phase (Fig. lc), is also seen to be true for orthorhombic and tetragonal 
phases: switching with fields applied along polar direction occurs at much higher fields 
than that applied along other crystallographic directions. While the polarization rotation {3} 
through mono clinic phases (which are known to be present near the morphotropic phase 



boundaries 



131 ]) has been shown to be responsible for giant piezoelectric response, 
we show that it facilitates polarization switching at lower fields and is crucial for memory 
applications. 

We further investigate the existence and thermodynamic stability of monoclinic phases 
through E — T phase diagrams. The electric field is turned on at the beginning of the sim- 
ulations, and is kept constant. The temperature is then increased (decreased) in order to 
perform heating (cooling) simulations which give the polarization, and dielectric suscepti- 
bility as a function of temperature. We simulated temperature dependent phase transitions 
for fields along [001], [101] and [111] directions varying from 10 kV/cm to 100 kV/cm in 
steps of 5 kV/cm. The contour plots of eigenvalues of dielectric tensor as a function of field 
and temperature naturally give a picture of the electric field - temperature (E — T) phase 
diagrams (Fig. 4). As the peak in dielectric constant marks a transition, high density of 
contours visible clearly as dark regions resembling thick lines/curves (Fig 4) mark the phase 
boundaries in E — T phase diagrams. These lines meet T-axis at 320 K, 250 K and 220 K at 
E = 0, which are the transitions temperatures in zero electric field [2 Oj. Hysteresis in heating 
and cooling simulations give small errors (~ 10 K and 10 kV/cm) in the phase boundaries . 
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FIG. 4: (Color online) E — T phase diagrams expressed as contour plots of eigenvalues of the 
dielectric tensor. High density of contours near the peak values of dielectric response give rise to 
thick and dark lines which form the boundaries between different stable phases. A crossover from 
one phase to another at such a boundary is also seen as a discontinuity in the slope of rj(E) in Fig. 
3. 

For field along [001] direction, there is no qualitative change in the phase transition 
behavior with increasing field (Fig 4a); the temperature range of the stable tetragonal phase 
expands with field. We note that a nonzero field breaks the symmetry of the system and 
the transition from C to T phase becomes diffused, and above the field of 50 kV/cm, there 
is virtually no difference between the T and C phases. The orthorhombic phase transforms 
to monoclinic Mq phase above a field of about 30 kV/cm, as characterized by a change 
in the slope of strain as a function of electric field (see Fig. 3). When the field is applied 
along the [101] direction, the temperature interval of the stable orthorhombic phase expands 
with field. The tetragonal phase transforms through an intermediate monoclinic Mc phase 
to an orthorhombic phase above 80v^2 kV/cm. Boundaries of the region of stability of 
the Mc phase are characterized by peaks in eigenvalues of dielectric response appearing as 
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thick dark lines (see Fig. 4b). Similarly for field along [111] direction, both tetragonal and 
orthorhombic phases transform at fields of about 40\/3 and 20a/3 kV/cm to Ma and Mb 
phases, respectively. Above the field of 60a/3 kV/cm, there is no clear boundary between 
rhombohedral and cubic phases. The lines of transition from monoclinic phases to T, O 
and R phases (with increasing fields) in our E — T phase diagrams are topologically similar 
to those in Bell's work[10] based on phenomenological theory and even the fields of these 
transitions agree within ~20% with Bell's estimates. However, our phase diagram differs 
in an important way: the region of stable monoclinic phases appears only above nonzero 
electric fields and its boundary with T, O and R phases at lower fields is characterized by 
(a) peak in the dielectric response (Fig. 4) and (b) change in the slope of rj changes (Fig. 
3). 

The sequence of phases appearing along the switching paths followed by the system during 
P-E hysteresis can be understood from the constant temperature lines in the E — T phase 
diagrams. For example, applying field along [111] direction of the tetragonal phase at 290 K 
would result in P tracing the path described in Fig. 2d. The connectivity between different 
monoclinic phases and T, O, R phases in the E — T phase diagram is consistent with that 
obtained from the Cohen's cube|3j. The change in the slope of strain as a function of field 
(see Fig 3) at a crossover from a T or O or R phase to a monoclinic one gives rise to the 
giant piezoelectric response arising at finite fields in experiments [7, 8]. 

To summarize, we have used effective hamiltonian in molecular dynamics simulations and 
have shown that: (1) switching fields are highest for the cases when the external electric 
field is applied along the direction of the polar axis, (2) switching occurs at low values of 
fields along other directions because of the polarization rotation through monoclinic phases. 
Stability of these monoclinic phases has been carefully analysed through determination of the 
E — T phase diagrams, which should be fundamental to the design of both electromechanical 
and memory type of devices. 
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